Analytic Approach to the Late Stages of Giant Planet Formation

The rate at which the circumplanetary disk receives infalling material from the envelope is given by the density field and the velocity field. Using the results from the previous section, we obtain

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d{\rm{\Sigma }}}{{dt}}=2\rho {v}_{\theta }\right|}_{(r,\mu =0)}=\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{2\pi {r}^{2}}\displaystyle \frac{{v}_{\theta }}{| {v}_{r}| }\displaystyle \frac{d{\mu }_{0}}{d\mu }=\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{4\pi {r}^{2}}\displaystyle \frac{1}{\zeta {\mu }_{0}},\end{eqnarray} \tag{ 20 }$$

where all quantities are evaluated at the disk plane (μ = 0). Note that the infall strikes the disk from both sides (top and bottom), leading to a factor of 2 in the result. Again using the orbit Equation (5), the above expression can be evaluated to obtain

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Sigma }}}{{dt}}=\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{4\pi {r}^{2}}\displaystyle \frac{1}{{\left({\zeta }^{2}-\zeta \right)}^{1/2}}.\end{eqnarray} \tag{ 21 }$$

As a benchmark, it is useful to determine the surface density of the disk that would result in the absence of disk accretion (see the following section). This surface density is given by the integral of dΣ/dt over time, equivalently mass, so that we obtain

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}=\int \displaystyle \frac{{\dot{M}}_{{\rm{p}}}{dt}}{4\pi {r}^{2}}\displaystyle \frac{1}{{\left({\zeta }^{2}-\zeta \right)}^{1/2}}.\end{eqnarray} \tag{ 22 }$$

The subscript indicates that this surface density would result in the absence of viscous evolution (ν = 0). Using ${dM}={\dot{M}}_{{\rm{p}}}{dt}$ and the definition of ζ,

$$\begin{eqnarray}&&\zeta \,=\,\displaystyle \frac{R}{3r}=\displaystyle \frac{a}{3r}{\left(\displaystyle \frac{M}{3{M}_{* }}\right)}^{1/3}\,\,\Rightarrow \,\,M={(3r/a)}^{3}3{M}_{* }{\zeta }^{3},\end{eqnarray} \tag{ 23 }$$

we find

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}=\displaystyle \frac{243{M}_{* }r}{4\pi {a}^{3}}\int \displaystyle \frac{{\zeta }^{2}d\zeta }{{\left({\zeta }^{2}-\zeta \right)}^{1/2}}.\end{eqnarray} \tag{ 24 }$$

Next we change variables to u defined such that u = 1/ζ, so that the limits of integration are u = 3r/R_H and u = 1. The resulting integral can be evaluated to obtain

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}=\displaystyle \frac{3{M}_{{\rm{p}}}}{16\pi {R}_{C}^{2}}\left\{\displaystyle \frac{\sqrt{1-u}(2+3u)}{u}+3u{\tanh }^{-1}(\sqrt{1-u})\right\}.\end{eqnarray} \tag{ 25 }$$

Once the circumplanetary disk has become well developed, the first term in the above expression dominates, and the surface density can be written in the approximate form

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}\approx \displaystyle \frac{3{M}_{{\rm{p}}}}{8\pi {R}_{C}r}.\end{eqnarray} \tag{ 26 }$$

The enclosed mass thus has the radial dependence M(r) ∝ r. As a result, most of the mass enters the planet–disk system at large radii r ≫ R_p, which in turn implies that disk accretion must be crucial for forming planets. ⁴

Figure 1 shows how the (projected) infalling trajectories of the flow join onto the circumplanetary disk. The inward paths are nearly radial at large distances r > R_C (note that the length scales in the figure are given in units of the centrifugal radius R_C = R_H /3). In the inner regime, however, the trajectories depart significantly from radial paths, crossing the midplane and joining the disk well before reaching the planetary surface. Most of the incoming trajectories thus intercept the disk rather than the planetary surface (which corresponds to x ≈ 0.006 in the figure).

For completeness, we note that the solution (25) for the surface density assumes that infall takes place over all initial polar angles θ₀. If infall is confined to the polar regions, as indicated by some numerical treatments (e.g., Lambrechts & Lega 2017; Lambrechts et al. 2019) then the centrifugal barrier R_C will be smaller and the resulting surface density distribution will be steeper (see Fung et al. 2019). This modification, considering only the range $0\leqslant {\theta }_{0}\leqslant {\theta }_{\max }$, is easily incorporated into this analytic treatment. On the other hand, viscous evolution will act to spread out the disk, increase its radius, and make the surface density distribution less steep, as shown in the following section.

Once the material has fallen onto the circumplanetary disk, its subsequent evolution occurs through viscous dissipation. Following standard nomenclature (Hartmann 2009), we write the viscosity in the form

$$\begin{eqnarray}&&\nu \,=\,\alpha {c}_{s}H=\alpha {c}_{s}^{2}/{\rm{\Omega }}.\end{eqnarray} \tag{ 27 }$$

Note that disk evolution for planet formation lies in a different regime than for star formation. The relevant quantity is the ratio of the viscous timescale to the accretion timescale, i.e.,

$$\begin{eqnarray}&&{ \mathcal R }\,=\,\displaystyle \frac{{R}^{2}}{\nu }\displaystyle \frac{\dot{M}}{M}\,=\,\displaystyle \frac{1}{\alpha }\displaystyle \frac{{R}^{2}}{{H}^{2}}\displaystyle \frac{\dot{M}P}{2\pi M}\sim \displaystyle \frac{100}{2\pi \alpha }\displaystyle \frac{\dot{M}P}{M},\end{eqnarray} \tag{ 28 }$$

where P is the orbit time at the outer disk edge and where the final approximate equality assumes H ∼ R/10. For star formation, R ∼ 100 au, M = M_* ∼ 1 M_⊙, P ∼ 1000 yr, and ${M}_{* }/\dot{M}\sim {10}^{5}$ yr. As a result, the ratio ${{ \mathcal R }}_{\mathrm{sf}}\sim 1/2\pi \alpha$. For planet formation, ${M}_{{\rm{p}}}/{\dot{M}}_{{\rm{p}}}\sim 1\,\mathrm{Myr}$, whereas P ∼ 1 yr, so that ${{ \mathcal R }}_{\mathrm{pf}}\sim {10}^{-4}/2\pi \alpha$. A much smaller viscosity is required for a planet-forming disk to keep up with the infall compared to a star-forming disk.

The equation of motion for a viscous accretion disk, including the infall terms derived in the previous section, can be written as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}=\displaystyle \frac{3}{r}\displaystyle \frac{\partial }{\partial r}\left[{r}^{1/2}\displaystyle \frac{\partial }{\partial r}\left({r}^{1/2}\nu {\rm{\Sigma }}\right)\right]+\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{4\pi {r}^{2}}\displaystyle \frac{1}{{\left({\zeta }^{2}-\zeta \right)}^{1/2}}.\end{eqnarray} \tag{ 29 }$$

In steady state, corresponding to vanishing time derivatives, the solution takes the form

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)=\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{6\pi \nu }\displaystyle \frac{1}{\sqrt{u}}\left[{\sin }^{-1}\sqrt{u}+\sqrt{u(1-u)}\right]=\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{6\pi \nu }f(u),\end{eqnarray} \tag{ 30 }$$

where u ≡ r/R_C and where the second equality defines the function f(u), which is of order unity and slowly varying with radius. Specifically, we note that the function f → 2 in the inner limit u → 0 and f → π/2 in the outer limit u → 1.

The diffusion equation has an additional solution, equivalently, an additional term, which has the form

$$\begin{eqnarray}&&\nu {\rm{\Sigma }}\,=\,\displaystyle \frac{{K}_{1}}{{r}^{1/2}}+{K}_{2},\end{eqnarray} \tag{ 31 }$$

where the K_j are integration constants, which can be chosen to satisfy the inner boundary condition. For stellar accretion disks, the usual assumption is to let νΣ → 0 at the inner boundary. If we use this inner boundary condition, the full solution can be written in the form

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)=\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{6\pi \nu }\left[f(u)-f({u}_{{\rm{p}}}){\left(\displaystyle \frac{{u}_{{\rm{p}}}}{u}\right)}^{1/2}\right],\end{eqnarray} \tag{ 32 }$$

where u_p ≡ R_p/R_C . The second term is negligible except near the surface of the planet. Moreover, in the present application, the disk is magnetically truncated so that the effective inner disk edge is well outside the planetary surface (see Section 4.3).

We can integrate the steady-state surface density over the disk area to find the disk mass:

$$\begin{eqnarray}&&{M}_{{\rm{d}}}=\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{3}{\int }_{0}^{{R}_{C}}\displaystyle \frac{{rdr}}{\nu }\displaystyle \frac{1}{\sqrt{u}}\left[{\sin }^{-1}\sqrt{u}+\sqrt{u(1-u)}\right].\end{eqnarray} \tag{ 33 }$$

If we assume that the temperature within the disk scales as T ∼ r^−1/2, then ν ∼ r, so that we can write ν = ν_C (r/R_C ) = ν_C u. The disk integral then becomes

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{\rm{d}}} & = & \displaystyle \frac{{\dot{M}}_{{\rm{p}}}{R}_{C}^{2}}{3{\nu }_{C}}{\displaystyle \int }_{0}^{1}\displaystyle \frac{{du}}{\sqrt{u}}\left[{\sin }^{-1}\sqrt{u}+\sqrt{u(1-u)}\right]\\ & = & \displaystyle \frac{{\dot{M}}_{{\rm{p}}}{R}_{C}^{2}}{3{\nu }_{C}}\left[\pi -\displaystyle \frac{4}{3}\right].\end{array}\end{eqnarray} \tag{ 34 }$$

We thus find that ${M}_{{\rm{d}}}/{M}_{{\rm{p}}}\sim { \mathcal R }$ from Equation (28), as expected.

With the steady-state solution in place, we can now consider how the disk evolves into such a configuration. The steady-state disk arises when the two terms on the right-hand side of Equation (29) are in balance. At the onset of disk formation, however, the viscous term is small (compared to the infall term) because the surface density must be small. In this regime, the surface density of the disk builds up mass from infall, as described in the previous section. The surface density will grow until it becomes comparable to the steady-state solution, where the required time for this transient growth is determined by the viscosity so that ${\rm{\Delta }}t\sim {R}_{C}^{2}/{\nu }_{C}$. Since the viscous timescale is much shorter than the evolutionary timescale (Equation (28)), the disk surface density has the steady-state form for most of its lifetime. For completeness, note that the centrifugal radius R_C grows with increasing planet mass with the scaling law ${R}_{C}\sim {M}_{{\rm{p}}}^{1/3}$. The weak dependence of R_C on mass implies that R_C varies only weakly with time, so that the evolution of disk surface density can be separated into the regimes described above.

With the outer boundary of the disk specified, we need to consider the inner boundary. In the absence of planetary magnetic fields, the circumplanetary disks can extend all the way to the planet's surface. In general, however, we expect the planet to develop a strong internal magnetic field with surface strength B_p (see Appendix B for an estimate). To leading order, such fields have a dipole form:

$$\begin{eqnarray}&&B(r,\theta )={B}_{{\rm{p}}}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{r}\right)}^{3}\left[3\cos \theta \hat{r}-\hat{z}\right],\end{eqnarray} \tag{ 35 }$$

where the spherical coordinate system is centered on the planet. The incoming flow from the circumplanetary disk cannot penetrate past the point where the ram pressure of the inward flow is balanced by the outward pressure from the magnetic field. This magnetic truncation radius has the form

$$\begin{eqnarray}&&{R}_{X}=\omega {\left(\displaystyle \frac{{B}_{{\rm{p}}}^{4}{R}_{{\rm{p}}}^{12}}{{{GM}}_{{\rm{p}}}{\dot{M}}_{{\rm{p}}}^{2}}\right)}^{1/7},\end{eqnarray} \tag{ 36 }$$

where ω is a dimensionless constant of order unity (Ghosh & Lamb 1978; Blandford & Payne 1982; see also the discussion of Mohanty & Shu 2008 for an alternate derivation).

If we insert typical numbers, the truncation radius can be evaluated to find

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{R}_{X}}{{R}_{{\rm{p}}}} & \approx & 3.8{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{J}}\right)}^{-1/7}{\left(\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{1\,{M}_{J}\,/\,\mathrm{Myr}}\right)}^{-2/7}\\ & & \times {\left(\displaystyle \frac{{B}_{{\rm{p}}}}{500\,\mathrm{gauss}}\right)}^{4/7}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{10}^{10}\,\mathrm{cm}}\right)}^{5/7}.\end{array}\end{eqnarray} \tag{ 37 }$$

For comparison, the centrifugal radius, and hence the disk radius, is given by

$$\begin{eqnarray}&&{R}_{{\rm{d}}}={R}_{C}=\displaystyle \frac{{R}_{H}}{3}\,\,\Rightarrow \,\,\displaystyle \frac{{R}_{{\rm{d}}}}{{R}_{{\rm{p}}}}\approx 170\ {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{J}}\right)}^{1/3}\left(\displaystyle \frac{a}{5\mathrm{au}}\right).\end{eqnarray} \tag{ 38 }$$

Since the disk radius scales as with the semimajor axis, planets forming in the inner regions of the circumstellar disk have radially small infall regions and produce radially small circumplanetary disks. For standard parameters, the magnetic truncation radius is larger than the disk radius for planets forming at semimajor axes a ≲ 0.1 au. As a result, magnetic fields could represent a significant obstacle to forming giant planets near their host stars. ⁵ For completeness, we can write the condition R_X = R_C in the form

$$\begin{eqnarray}&&\left(\displaystyle \frac{a}{{R}_{{\rm{p}}}}\right){\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{* }}\right)}^{10/21}\,=\,{3}^{4/3}\omega {\left(\displaystyle \frac{{B}_{{\rm{p}}}^{4}{R}_{{\rm{p}}}^{5}}{{{GM}}_{* }{\dot{M}}_{{\rm{p}}}^{2}}\right)}^{1/7}.\end{eqnarray} \tag{ 39 }$$

The previous subsections defined the expected structure for the circumplanetary disks that arise during the process of planet formation. The surface densities for these disks are shown in Figure 2. The upper black dashed curve shows the surface density profile that would result in the absence of viscous accretion (as derived in Section 4.1). The surface density shows a radial dependence of the approximate form Σ ∝ 1/r, except near the outer edge, where an enhancement is present, followed by a sharp edge at r = R_C . In the presence of viscosity, the disk surface density is given by the steady-state solutions of Section 4.2. Surface density profiles are shown in the figure for α = 10⁻², 10⁻³, and 10⁻⁴, from bottom to top. The larger viscosity results in lower disk masses in steady state. The surface density profiles are normalized by dividing out a factor of ${M}_{{\rm{p}}}/{R}_{C}^{2}$, where M_p represents the total mass that has fallen to the planet–disk system. For the viscous models, we further specify the mass infall rate to be 1 M_J /Myr and the radial location of the planet to be a = 5 au. With these choices, for a planet with mass M_p = 1 M_J , the circumplanetary disk has outer radius R_C /R_p = 170 and an inner boundary given by the magnetic truncation radius at R_X /R_p = 3.8, as depicted by the vertical orange dashed line in Figure 2. Although the surface density profiles are truncated at r = R_C in the figure, the solutions can be extended beyond this nominal disk radius where they match onto more steeply decreasing profiles (see Appendix D). Note that the disk must spread beyond R_C in order to satisfy conservation of angular momentum.

Zoom In Zoom Out Reset image size

The circumplanetary disks thus have relatively simple forms. The dynamic range is small, with R_C /R_X ∼ 40−50. The disk masses are a small fraction of the total mass (Equation (34)) and the surface density distribution is close to a 1/r power law (Equation (30)). In spite of the small disk masses, the surface densities are large enough for the disks to remain optically thick to their internal infrared radiation.

Before leaving this section, it is useful to consider how gaseous planets would form in the absence of disk accretion (through the circumplanetary disk). If the viscosity vanishes, ν → 0, then the rate at which the planet gains mass is limited by the geometry of the infall and has the form

$$\begin{eqnarray}&&{\left.{\dot{M}}_{{\rm{p}}}\right|}_{\mathrm{direct}}\,=\,4\pi {R}_{{\rm{p}}}^{2}{\int }_{0}^{1}\rho | {v}_{r}| d\mu ,\end{eqnarray} \tag{ 40 }$$

where the integrand is evaluated at the planetary surface. Using the infall solution from Section 3 (see also Adams & Shu 1986, for the star formation analog), we can evaluate the integral to obtain

$$\begin{eqnarray}&&{\left.{\dot{M}}_{{\rm{p}}}\right|}_{\mathrm{direct}}\,=\,{\dot{M}}_{{\rm{p}}}\left[1-\sqrt{1-{R}_{{\rm{p}}}/{R}_{C}}\right].\end{eqnarray} \tag{ 41 }$$

The rate of accretion onto the planet is thus an ever-decreasing fraction of the total infall rate.

The planet mass, M_p, for a given mass, M, that has fallen is given by the integral expression

$$\begin{eqnarray}&&{M}_{{\rm{p}}}={M}_{0}+{\int }_{{M}_{0}}^{M}{dM}\left[1-\sqrt{1-{({M}_{0}/M)}^{1/3}}\right],\end{eqnarray} \tag{ 42 }$$

where the scale M₀ is the mass for which the centrifugal radius is equal to the planetary radius, i.e.,

$$\begin{eqnarray}&&{M}_{0}=81{M}_{* }{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{a}\right)}^{3}.\end{eqnarray} \tag{ 43 }$$

For planets forming at a ∼ 5 au, M₀ ∼ 2 × 10⁻⁷ M_J . In contrast, for hot Jupiters at a = 0.05 au, M₀ ∼ 0.2 M_J . For both scenarios, the formation of a Jovian-mass planet requires accretion via a circumplanetary disk. For cold Jupiters, the vast majority of the mass must be accreted through the disk. For hot Jupiters, however, a significant fraction of the mass can be accreted directly (notice also that hot Jupiters tend to have somewhat smaller masses).

Note that the above treatment does not include the effects of planetary magnetic fields, which act to increase the effective-capture cross section of the growing planet. For a given estimate of the effective-capture radius R_σ (see Section 4.6) one can replace R_p with R_σ in the above expressions.

A related question is: How much total mass must flow through the Hill sphere in order for a given mass to directly accrete onto the planet? We can evaluate Equation (42) in the asymptotic limit M ≫ M₀ to obtain

$$\begin{eqnarray}&&{M}_{{\rm{p}}}\to \displaystyle \frac{3}{4}{\left({M}_{0}{M}^{2}\right)}^{1/3}.\end{eqnarray} \tag{ 44 }$$

In order to produce a final planet mass of M_p = 1 M_J (with no circumplanetary disk accretion and at radial location a = 5 au), we would need a total mass M ≈ 3400 M_J = 3.4 M_⊙. In other words, in order for a Jovian mass to fall directly to the planet, several stellar masses of material must cycle through the Hill sphere. Some mechanism for redistributing angular momentum—within the sphere of influence of the forming object—is thus necessary for the formation of Jovian planets. Viscosity in the circumplanetary disk provides such a mechanism for the scenario explored in this paper. ⁶

The magnetic truncation radius R_X (Equation (36)), which arises from the balance between disk accretion and magnetic pressure, defines the inner edge of the circumplanetary disk. However, for material falling onto the system near the rotational poles, a second magnetic boundary, R_σ , arises from the balance between the ram pressure of infalling material and the magnetic pressure. These two scales are different because not all of the infalling material can reach inner radii comparable to the planetary size, so that the mass infall rate is reduced and the boundary is larger. Equating the incoming ram pressure ρ v² with the magnetic field pressure at the pole (μ = 1), we can write the magnetic capture radius in the form

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{\sigma }}{{R}_{{\rm{p}}}}\,=\,\tilde{\omega }{\left(\displaystyle \frac{{B}_{{\rm{p}}}^{4}{R}_{{\rm{p}}}^{3}{R}_{C}^{2}}{{{GM}}_{{\rm{p}}}{\dot{M}}_{{\rm{p}}}^{2}}\right)}^{1/9},\end{eqnarray} \tag{ 45 }$$

where $\tilde{\omega }$ is a dimensionless constant of order unity and where the quantity ${\dot{M}}_{{\rm{p}}}$ is the total mass accretion rate. For typical parameter values, the capture radius is somewhat larger than the disk truncation radius so that R_σ ∼ 1.4 R_X (compare with Equation (36)).

The magnetic field lines and incoming trajectories near the forming planet are shown in Figure 3. Keep in mind that all of the trajectories in the figure correspond to small values of the initial polar angle θ (otherwise the trajectories would cross the plane of the disk far to the right). In order for an incoming parcel of gas to be deflected by the magnetic fields, its radial location is given by Equation (45). In order for the parcel to subsequently move inward onto the planet, rather than outward and back toward the disk, it must fall sufficiently close to the pole. Here we approximate this bifurcation point as the point where the magnetic field line is horizontal and where the magnetic pressure balances the ram pressure. This critical point of the field line corresponds to the location where the horizontal coordinate has the form

$$\begin{eqnarray}&&{x}_{\mathrm{hor}}\,=\,{\left(\displaystyle \frac{2}{3}\right)}^{1/2}{R}_{\sigma }\sim 0.82{R}_{\sigma }.\end{eqnarray} \tag{ 46 }$$

We thus have the following result: the capture radius R_σ is larger than the magnetic truncation radius because the direct infall is reduced in the inner region. However, only the innermost portion of the capture region funnels material directly onto the planet rather than outward and onto the disk (see Batygin 2018 and Batygin & Morbidelli 2020). The radius R_σ is somewhat larger than R_X and x_hor is somewhat smaller than R_σ , so that the two corrections tend to compensate. In the end, the effective-capture cross section is given approximately by the length scale R_X .

Zoom In Zoom Out Reset image size

Since most of the mass falling onto the planet–disk system is eventually accreted onto the growing planet, the luminosity contributions can be written in terms of the benchmark power scale

$$\begin{eqnarray}&&{L}_{0}\equiv \displaystyle \frac{{{GM}}_{{\rm{p}}}{\dot{M}}_{{\rm{p}}}}{{R}_{{\rm{p}}}},\end{eqnarray} \tag{ 47 }$$

which corresponds to the total power that could be generated if all of the incoming material falls (from rest at infinity) onto the planetary surface at freefall speeds and dissipates all of its kinetic energy.

The luminosity of the system is distributed between the planet and its disk. Only a fraction, f_p, of the incoming material falls directly onto the planet (see Equation (41)). However, an additional fraction of the material falls within the magnetic capture radius (Section 4.6) and is funneled onto the planet along magnetic field lines. The remaining incoming material falls onto the disk. For the sake of definiteness and simplicity, we assume that the magnetic truncation radius R_X of the disk and the magnetic capture radius R_σ are equal. The fractions of the incoming material that fall onto the planet and onto the disk are thus given by

$$\begin{eqnarray}&&{f}_{{\rm{p}}}\,=\,1-\sqrt{1-{u}_{{\rm{x}}}}\quad \mathrm{and}\quad {f}_{{\rm{d}}}\,=\,\sqrt{1-{u}_{{\rm{x}}}},\end{eqnarray} \tag{ 48 }$$

where we have defined u_x ≡ R_X /R_C . If the planetary magnetic field is sufficiently weak, so that R_X ≤ R_p, then u_x = R_p/R_C .

To specify the planetary component of the accretion luminosity, we assume that all of the material that falls within the capture radius is accreted by the planet, and that shocks on the surface dissipate the radial and poloidal components of the velocity. The planet is assumed to be corotating with its disk at the magnetic truncation point, so that the planetary rotation rate is given by

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{p}}}^{2}=\displaystyle \frac{{{GM}}_{{\rm{p}}}}{{R}_{X}^{3}}.\end{eqnarray} \tag{ 49 }$$

The incoming material thus retains a small azimuthal velocity given by ${v}_{\phi }\,=\,\sin \theta {R}_{{\rm{p}}}{{\rm{\Omega }}}_{{\rm{p}}}$. The luminosity thus takes the form

$$\begin{eqnarray}&&{L}_{{\rm{p}}}^{\mathrm{dir}}\,=\,\left(1-\sqrt{1-{u}_{{\rm{x}}}}\right)\left(1-\displaystyle \frac{1}{2}\langle {\sin }^{2}\theta \rangle \displaystyle \frac{{R}_{{\rm{p}}}^{3}}{{R}_{X}^{3}}\right){L}_{0},\end{eqnarray} \tag{ 50 }$$

where the angular brackets represent an average over the colatitudes. The newly added material is eventually distributed over the planetary surface, so we expect $\langle {\sin }^{2}\theta \rangle \approx 2/3$.

Although the remaining fraction of the material falls onto the disk, it eventually accretes onto the planet. However, it falls from the magnetic truncation radius, rather than from infinity, so that the additional planetary luminosity derived from accretion takes the approximate form

$$\begin{eqnarray}&&{L}_{{\rm{p}}}^{\mathrm{acc}}\,=\,{L}_{0}\sqrt{1-{u}_{{\rm{x}}}}\left(1-\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{X}}\right)\left(1-\displaystyle \frac{{R}_{{\rm{p}}}^{3}}{3{R}_{X}^{3}}\right).\end{eqnarray} \tag{ 51 }$$

As written, the expression has three correction factors. The first determines the amount of material that falls onto the disk (the remaining fraction is accounted for in the previous contribution). The second factor takes into account the starting point for the material, which falls from an initial radius R_X instead of r → ∞ . The third factor takes into account the rotational motion of the incoming material, where we again assume that the magnetic field causes the incoming material to rotate with the mean motion of the disk at the location R_X and that the incoming material is uniformly distributed in planetary latitude. Although this correction factor is close to unity, note that other assumptions can be invoked. For completeness, we note that the planet and its disk could support outflows or winds analogous to those observed in young stellar objects. Such winds could result in a reduced accretion rate onto the planetary surface, and could play a role in redistributing angular momentum. Although potentially important, we leave this issue for future work.

The luminosity contributions from both direct infall and from the inner edge of the disk combine to determine the total planetary accretion luminosity L_p, which takes the form

$$\begin{eqnarray}&&{L}_{{\rm{p}}}\,=\,\left(1-\sqrt{1-{u}_{{\rm{x}}}}\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{X}}\right)\left(1-\displaystyle \frac{{R}_{{\rm{p}}}^{3}}{3{R}_{X}^{3}}\right){L}_{0}.\end{eqnarray} \tag{ 52 }$$

In general, the planet has an additional internal luminosity, L_int, although the accretion luminosity is expected to be somewhat larger (see, e.g., Marley et al. 2007). The planet must radiate the total luminosity, given by Equation (52) and L_int, and has an effective surface temperature given by

$$\begin{eqnarray}&&\sigma {T}_{{\rm{p}}}^{4}=\displaystyle \frac{{L}_{0}}{4\pi {R}_{{\rm{p}}}^{2}}\left(1-\sqrt{1-{u}_{{\rm{x}}}}\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{X}}\right)\left(1-\displaystyle \frac{{R}_{{\rm{p}}}^{3}}{3{R}_{X}^{3}}\right)+\displaystyle \frac{{L}_{\mathrm{int}}}{4\pi {R}_{{\rm{p}}}^{2}}.\end{eqnarray} \tag{ 53 }$$

In the limit L_int → 0, only the first term contributes, and we make this approximation for the rest of this paper.

The remainder of the energy is dissipated within the disk, which has total luminosity given by

$$\begin{eqnarray}&&{L}_{{\rm{d}}}\,=\,{L}_{0}\displaystyle \frac{{R}_{{\rm{p}}}}{2{R}_{X}}\sqrt{1-{u}_{{\rm{x}}}},\end{eqnarray} \tag{ 54 }$$

where the factor of 2 arises because half of the energy is stored in the rotational energy of the disk. Because the magnetic truncation radius is typically R_X ∼ 4R_p, the total accretion luminosity generated on the planetary surface is larger than the disk luminosity by (about) an order of magnitude. From the perspective of the infalling envelope, most of the energy is generated by a point source in the center of the structure.

The total disk luminosity is distributed among three components. The perpendicular component, v_θ , of the incoming velocity is dissipated as it strikes the disk, thereby producing luminosity ${L}_{{\rm{d}}}^{\mathrm{surf}}$. When the material reaches the disk, however, it generally will not have the right azimuthal velocity to enter into Keplerian rotation about the planet. The new material thus mixes with existing disk material, dissipating energy while conserving angular momentum, and produces a mixing luminosity, ${L}_{{\rm{d}}}^{\mathrm{mix}}$. The remaining luminosity, ${L}_{{\rm{d}}}^{\mathrm{visc}}$, must be dissipated (in this approximation) by the viscosity. The surface and mixing contributions for circumstellar disks have been calculated previously (Cassen & Moosman 1981; Cassen & Summers 1983; Adams & Shu 1986) and can be written in the form

$$\begin{eqnarray}\begin{array}{rcl}{L}_{{\rm{d}}}^{\mathrm{surf}} & = & {\displaystyle \int }_{A}\left(\displaystyle \frac{1}{2}{v}_{\theta }^{2}\right)\rho | {v}_{\theta }| {dA}=\displaystyle \frac{{u}_{{\rm{x}}}}{4}\\ & & \times \left\{\mathrm{ln}\left[\displaystyle \frac{1+\sqrt{1-{u}_{{\rm{x}}}}}{1-\sqrt{1-{u}_{{\rm{x}}}}}\right]-2\sqrt{1-{u}_{{\rm{x}}}}\right\}{L}_{0},\end{array}\end{eqnarray} \tag{ 55 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{L}_{{\rm{d}}}^{\mathrm{mix}} & = & {\displaystyle \int }_{A}\displaystyle \frac{1}{2}\left\{{v}_{r}^{2}+{\left[{v}_{\phi }-{\left(\displaystyle \frac{{{GM}}_{{\rm{p}}}}{r}\right)}^{1/2}\right]}^{2}\right\}\rho | {v}_{\theta }| {dA}\\ & = & \displaystyle \frac{{u}_{{\rm{x}}}}{2}\left\{\mathrm{ln}\left[\displaystyle \frac{1+\sqrt{1-{u}_{{\rm{x}}}}}{1-\sqrt{1-{u}_{{\rm{x}}}}}\right]+\sqrt{1-{u}_{{\rm{x}}}}+\pi \right.\\ & & -\left.2{\tan }^{-1}{\left[\displaystyle \frac{{u}_{{\rm{x}}}}{1-{u}_{{\rm{x}}}}\right]}^{1/2}\right\}{L}_{0}.\end{array}\end{eqnarray} \tag{ 56 }$$

These two luminosity contributions (see also Szulágyi & Mordasini 2017) represent the energy that is dissipated as material joins the disk, independent of any subsequent viscous evolution, and can be combined to define a total direct luminosity contribution ${L}_{{\rm{d}}}^{\mathrm{dir}}$ = ${L}_{{\rm{d}}}^{\mathrm{surf}}+$ ${L}_{{\rm{d}}}^{\mathrm{mix}}$, i.e.,

$$\begin{eqnarray}&&{L}_{{\rm{d}}}^{\mathrm{dir}}\,=\,\displaystyle \frac{{u}_{{\rm{x}}}}{4}\left\{3\mathrm{ln}\left[\displaystyle \frac{1+\sqrt{1-{u}_{{\rm{x}}}}}{1-\sqrt{1-{u}_{{\rm{x}}}}}\right]+2\pi -4{\sin }^{-1}\sqrt{{u}_{{\rm{x}}}}\right\}{L}_{0}.\end{eqnarray} \tag{ 57 }$$

In the limit where all of the material can be accreted onto the planet, the remaining viscous luminosity is given by

$$\begin{eqnarray}&&{L}_{{\rm{d}}}^{\mathrm{visc}}\,=\,{L}_{{\rm{d}}}-{L}_{{\rm{d}}}^{\mathrm{dir}}.\end{eqnarray} \tag{ 58 }$$

Although we can separate the contributions to the disk luminosity, the net result of viscous evolution is to transport all of the incoming material inward to the magnetic truncation radius. As a result, the total disk luminosity has the simple form given by Equation (54). We can approximate the temperature distribution from the disk surface using the form

$$\begin{eqnarray}&&{T}_{{\rm{d}}}(r)={T}_{X}{\left(\displaystyle \frac{{R}_{X}}{r}\right)}^{3/4},\end{eqnarray} \tag{ 59 }$$

where the leading coefficient is determined by the requirement that the disk radiates its total luminosity:

$$\begin{eqnarray}&&\sigma {T}_{X}^{4}=\displaystyle \frac{{{GM}}_{{\rm{p}}}{\dot{M}}_{{\rm{p}}}}{8\pi {R}_{X}^{3}}\displaystyle \frac{1}{\sqrt{1-{u}_{{\rm{x}}}}}.\end{eqnarray} \tag{ 60 }$$

Note that the power-law dependence of the temperature distribution (59) does not continue to arbitrarily large radii. The forming planet is embedded within the circumstellar disk of its host star, and this nebula provides a background minimum temperature for the circumplanetary disk.

With the surface temperatures of the planet and the disk specified, we can define the corresponding spectral energy distributions of the two central components:

$$\begin{eqnarray}&&{L}_{{\rm{p}}\nu }\,=\,4\pi {R}_{{\rm{p}}}^{2}\pi {B}_{\nu }({T}_{{\rm{p}}}),\end{eqnarray} \tag{ 61 }$$

and

$$\begin{eqnarray}&&{L}_{{\rm{d}}\nu }\,=\,2(2\cos \theta ){\int }_{{R}_{X}}^{{R}_{C}}\pi {B}_{\nu }[{T}_{{\rm{d}}}(r)]2\pi {rdr}.\end{eqnarray} \tag{ 62 }$$

The leading factor of 2 arises because the disk has two sides and the factor of $(2\cos \theta )$ arises from the viewing angle. Here we have assumed that the disk is optically thick to its own emission and this expression does not include external extinction (see also Zhu 2015 and Zhu et al. 2018). The corresponding monochromatic fluxes arising from the planet and its disk have the form

$$\begin{eqnarray}&&{F}_{{\rm{p}}\nu }(r,\mu )=\displaystyle \frac{{L}_{{\rm{p}}\nu }}{4\pi {r}^{2}}{G}_{\mathrm{pd}}(\mu )\exp [-{\tau }_{\nu }(r,\mu )],\end{eqnarray} \tag{ 63 }$$

and

$$\begin{eqnarray}&&{F}_{{\rm{d}}\nu }(r,\mu )=\displaystyle \frac{{L}_{{\rm{d}}\nu }}{4\pi {r}^{2}}(2\mu ){G}_{\mathrm{dp}}(\mu )\exp [-{\tau }_{\nu }(r,\mu )],\end{eqnarray} \tag{ 64 }$$

as seen by an observer with location (r, μ). The optical depth of the infalling envelope, τ_ν , depends on the viewing angle. For completeness, we include the functions G_dp(μ) and G_pd(μ), which take into account the mutual shadowing of the planet by the disk and the disk by the planet (Adams & Shu 1986). Most of the shadowing takes place near the planet, however, and we expect the magnetic truncation radius to remove the nearest disk material, so that the shadowing factors are close to unity in the present application.

Here we consider the infalling envelope surrounding the forming planet to be optically thin to its own emitted thermal photons. In this limit, the spectral energy, ${{ \mathcal L }}_{\nu }$, corresponding to dust emission has the form

$$\begin{eqnarray}&&{{ \mathcal L }}_{\nu }\,=\,{f}_{e}{\int }_{{R}_{X}}^{{R}_{H}}\rho {\kappa }_{\nu }4\pi {B}_{\nu }(T)4\pi {r}^{2}{dr},\end{eqnarray} \tag{ 65 }$$

where κ_ν is the frequency-dependent opacity and B_ν is the Planck function. The magnetic truncation radius R_X delineates the inner boundary of the envelope and the Hill radius R_H marks the outer boundary of the infall region. We are thus implicitly assuming that the radiation field reaches its eventual free-streaming form for r ≲ R_H . The factor f_e ≤ 1 represents a covering fraction, which allows for a possible reduction in the solid angle subtended by the envelope, with a corresponding reduction in the emitted radiation (although we take f_e = 1 for the sake of definiteness). The density distribution is that found in Section 3, where we use the equivalent spherical distribution for simplicity. ⁷

The interaction between the radiation field and the envelope is dominated by the dust opacity κ_ν . Although the dust opacity in the interstellar medium has been well characterized (starting with Draine & Lee 1984), the dust within circumstellar disks is expected to evolve significantly from its interstellar form. Dust grains generally grow, and opacity decreases as the grain size becomes comparable to the wavelength of the radiation. In the present context, the temperature of the envelope falls in the range 100 K ≤ T ≤ 1000 K. Over this interval, the Planck mean and Rosseland mean opacities are slowly varying functions of temperature, and the different possible models are in relatively good agreement (e.g., see Semenov et al. 2003, and references therein). As a result, to a good approximation, we can consider the opacity to have the power-law form ${\kappa }_{\nu }\,=\,{\kappa }_{0}{(\nu /{\nu }_{0})}^{\eta }$. The Planck mean opacity then has the temperature dependence κ_P = b_κ T^η , where the constant b_κ is given by

$$\begin{eqnarray}&&{b}_{\kappa }\equiv {\kappa }_{0}{\left(\displaystyle \frac{k}{h{\nu }_{0}}\right)}^{\eta }\displaystyle \frac{{\rm{\Gamma }}(4+\eta ){\zeta }_{R}(4+\eta ))}{6\,{\zeta }_{R}(4)}.\end{eqnarray} \tag{ 66 }$$

Here k and h are the Boltzmann and Planck constants, respectively, Γ(z) is the gamma function and ζ_R (z) is the Riemann zeta function (Abramowitz & Stegun 1972). With the choices η = 1, κ₀ = 10 cm²/g, and ν₀ = 10¹⁴ Hertz, the resulting κ_ν is somewhat smaller than the interstellar opacity, but the resulting mean opacities are within the ranges advocated by Semenov et al. (2003) for circumstellar disks. We use these values for the results presented below, but they are readily varied given the analytical solutions developed herein.

As a further simplification, we note that the infalling envelope is optically thin to its emitted radiation field. In this regime, we can solve the radiative transfer equation (Adams & Shu 1985) and find that the temperature profile takes the power-law form

$$\begin{eqnarray}&&T(r)={T}_{C}{\left(\displaystyle \frac{r}{{R}_{C}}\right)}^{-2/5},\end{eqnarray} \tag{ 67 }$$

where the power-law exponent is appropriate for the dust opacity law (κ_ν ∝ ν) used above (in general, the exponent is 2/(η + 4)). With the temperature profile specified, we can evaluate the total luminosity emitted by the infalling envelope to find

$$\begin{eqnarray}&&{ \mathcal L }\,=\,16\pi {\int }_{{R}_{X}}^{{R}_{H}}\rho {\kappa }_{P}\sigma {T}^{4}{r}^{2}{dr},\end{eqnarray} \tag{ 68 }$$

where κ_P is the Planck mean opacity. The integral expression for the envelope luminosity, ${ \mathcal L }$, can be evaluated and has the form

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal L } & = & 16\pi \,\sigma \,{{Cb}}_{\kappa }{T}_{C}^{5}{R}_{C}^{3/2}{\int }_{{u}_{{\rm{x}}}}^{3}{u}^{-3/2}{ \mathcal A }(u){du}\\ & = & 16\pi {R}_{C}^{2}\,\sigma \,{b}_{\kappa }{T}_{C}^{5}{N}_{\mathrm{col}}.\end{array}\end{eqnarray} \tag{ 69 }$$

Note that the upper limit of the integral in Equation (69) is taken to be the outer boundary of the infalling envelope at u = 3 (equivalently, r = R_H ). ⁸ The fiducial temperature, T_C , is thus given by

$$\begin{eqnarray}&&{T}_{C}={\left(\displaystyle \frac{{ \mathcal L }}{16\pi {R}_{C}^{2}\sigma {N}_{\mathrm{col}}{b}_{\kappa }}\right)}^{1/5}.\end{eqnarray} \tag{ 70 }$$

The total envelope luminosity ${ \mathcal L }$ is determined by the fraction of the central source radiation that is attenuated by the infalling envelope,

$$\begin{eqnarray}&&{ \mathcal L }\,=\,\left[{L}_{{\rm{p}}}+{L}_{{\rm{d}}}-{\int }_{0}^{\infty }d\nu ({L}_{{\rm{p}}\nu }+{L}_{{\rm{d}}\nu })\exp (-{\tau }_{\nu })\right],\end{eqnarray} \tag{ 71 }$$

where τ_ν = κ_ν N_col is the total optical depth through the envelope for a given frequency. For typical parameters, the fraction of the central source luminosity that is absorbed by the envelope is ∼ 0.25 − 0.5, i.e., ${ \mathcal L }\sim 0.25$ (L_p + L_d). The optical depth of the envelope to the central source radiation is thus of order unity. The optical depth of the envelope to its emitted radiation, which has longer wavelength, is much less than unity, so that the optically thin approximation is valid.

With the properties of the planet, disk, and infalling envelope specified, we can determine the spectral energy distribution of the forming planet. Here we present only a simple continuum model where the planetary surface and the disk surface are assumed to radiate as blackbodies. Future work should consider line emission.

The spectral energy distributions of the central planet–disk system are shown in Figure 4. The planet mass is taken to be M_p = 1 M_J with a mass accretion rate ${\dot{M}}_{{\rm{p}}}$ = 1 M_J /Myr. The planet is assumed to form at radial location a = 5 au, which determines its Hill radius and hence the disk radius. For a given choice of magnetic field strength on the planetary surface, the inner edge of the disk is determined by the magnetic truncation radius. The curves shown in the figure correspond to logarithmically spaced values in the range B_p = 10–1000 Gauss (increasing from bottom to top). Stronger magnetic fields result in larger magnetic truncation radii R_X . Under the assumption of steady-state disk accretion, larger values of R_X result in larger overall luminosity and more accretion power generated on the planetary surface (see Section 5.1). Smaller inner radii R_X allow for relatively more energy to be dissipated in the disk, rather than on the planetary surface, so that the spectral energy distributions become somewhat redder. The differences are modest, however, with the total luminosity varying by only a factor of ∼2 over the given parameter range.

Zoom In Zoom Out Reset image size

Radiation from the central planet–disk system is absorbed and reradiated by the infalling envelope (Section 5.2). The total optical depth of the envelope to the central source photons is less than unity, so that only a fraction of the total luminosity is reprocessed. Figure 5 shows the resulting emergent spectral energy distribution for a typical system with mass M_p = 1 M_J , mass accretion rate $\dot{M}$ = 1 M_J /Myr, location a = 5 au, and magnetic field strength B_p = 1000 Gauss (so that R_X ≈ 5.6 R_p). The contributions from the central source and the envelope are shown separately (dashed curves), with the total given by the solid curve. As expected, the envelope dominates the emission at sufficiently long wavelengths (λ ∼ 10 μm), whereas more energy is emitted at shorter wavelengths (λ ≲ 1 μm).

Zoom In Zoom Out Reset image size

The emergent spectral energy distributions depend on four variables that define the system: the planet mass M_p, the mass accretion rate ${\dot{M}}_{{\rm{p}}}$, the magnetic field strength B_p on the planetary surface, and the semimajor axis a of the planetary orbit. With $({M}_{{\rm{p}}},{\dot{M}}_{{\rm{p}}})$ specified, the field strength determines the magnetic truncation radius, which sets the inner edge of the disk. The semimajor axis determines the centrifugal radius, which sets the outer edge of the disk. We note that if a complete theoretical description of the planet formation process were available, then the dependence of the mass accretion rate and the magnetic field strength on planet mass and location could be calculated (see also Appendix E). ⁹ In the absence of a complete theory, it is useful to see how the spectral energy distributions vary over the parameter space $({M}_{{\rm{p}}},{\dot{M}}_{{\rm{p}}},{B}_{{\rm{p}}},a)$. For the sake of definiteness, we fix the planet radius so that R_p = 10¹⁰ cm.

Figure 6 shows how the spectral energy distributions depend on the planet mass, which is varied over the range M_p = 0.1–10 M_J . For all of the models, the mass accretion rate is fixed at ${\dot{M}}_{{\rm{p}}}$ = 1 M_J /Myr, the semimajor axis a = 5 au, and surface magnetic field strength B_p = 500 Gauss. Two trends are evident. First, as expected, the total luminosity of the planet–disk system increases with mass, so that the larger objects are brighter. Moreover, to leading order, L ∝ M_p. Second, the spectral energy distributions become bluer with increasing planet mass. For planets of higher mass, the planetary surface is hotter and the total optical depth through the infalling envelope is lower, so that a smaller fraction of the energy is emitted at longer wavelengths. These trends are consistent with previous numerical studies (Szulágyi et al. 2019).

Zoom In Zoom Out Reset image size

The dependence of the spectral energy distributions on the mass accretion rate is shown in Figure 7, where ${\dot{M}}_{{\rm{p}}}$ = 0.1–10 M_J /Myr. In this case, the planet mass is fixed at M_p = 1 M_J , the semimajor axis a = 5 au, and surface magnetic field strength B_p = 500 Gauss. As before, two trends are present (compare with Figure 6). The total luminosity of the planet–disk system increases with the mass accretion rate, where the luminosity scales as $L\propto {\dot{M}}_{{\rm{p}}}$ to leading order. In this case, however, the spectral energy distributions become redder with increasing luminosity (due to the larger mass accretion rate). The increase in ${\dot{M}}_{{\rm{p}}}$ corresponds to an increase in the column density, so that more of the central source radiation is absorbed by the infalling envelope and reradiated into the infrared.

Zoom In Zoom Out Reset image size

The effects of varying the magnetic field strength are shown in Figure 8. In this case, the B_p varies over the range 10–1000 Gauss. The planet has mass M_p = 1 M_J , accretion rate ${\dot{M}}_{{\rm{p}}}=1{M}_{J}$/Myr, and semimajor axis a = 5 au. For the smallest field value, B_p = 10 Gauss, the magnetic truncation radius R_X is smaller than the planetary radius, so that the disk extends all the way to the planetary surface. For stronger magnetic fields with B_p = 100 and 1000 Gauss, the radius R_X /R_p = 1.5 and 5.6, respectively, and magnetic truncation becomes important. The larger values of R_X lead to modest increases in the luminosity, due to less rotational energy being stored in the circumplanetary disk. Larger truncation radii R_X also result in bluer planet–disk spectra, so that the extinction (for fixed column density) is greater. As a result, a double-humped form of the spectral energy distribution is produced for systems with large field strengths B_p.

Zoom In Zoom Out Reset image size

Figure 9 shows how the spectral energy distributions vary with the semimajor axis of the orbit of the forming planet. Here the planet mass M_p = 1 M_J , the accretion rate ${\dot{M}}_{{\rm{p}}}$ = 1 M_J /Myr, and the magnetic field strength is fixed at B_p = 350 Gauss. Results are shown for semimajor axes in the range 0.1–10 au. The spectral energy distributions become wider as the semimajor axis increases. The most important effect is that the radius of the circumplanetary disk scales as R_d ∝ R_H ∝ a. For the lower end of the orbital range (a = 0.1 au), the disk radius becomes comparable to the magnetic truncation radius. In this regime, most of the radiation from the central source is emitted from a single surface temperature. Moreover, the column density is larger, so that more of the central source photons are absorbed and reemitted by the envelope. The resulting spectral energy distribution thus approaches a blackbody form. For larger semimajor axes, the column density is smaller, and the spectral energy distributions are broader and bluer.

Zoom In Zoom Out Reset image size

Our main results can be summarized as follows.

We have developed an infall solution to describe the density and velocity fields of the material falling toward a forming giant planet. The flow enters into the sphere of influence of the planet near the Hill radius, and then smoothly approaches ballistic trajectories. In this approximation, the density distribution ρ(r, μ) (given by Equations (5) and (12)) and the velocity field v (r, μ) (given by Equations (8)–(11)) are axisymmetric. The material with the highest specific angular momentum falls onto a circumplanetary disk with radius given by the centrifugal barrier (R_d = R_C = R_H /3; Equation (7)).

To gain insight into the nature of the region surrounding the forming planet, we have constructed the equivalent spherical density distribution (Equations (13)–(15)). The total column density N_col (given by Equation (17)) of the envelope plays an important role in determining the spectral energy distribution of the forming body. The envelope surrounding the planet is predicted to be marginally optically thick to the radiation emitted from the central planet–disk system, but is predicted to be optically thin to radiation emitted by the envelope itself.

Most of the material falling toward the planet has too much specific angular momentum to reach the planetary surface. Instead, material falls onto a circumplanetary disk with radius R_C determined by conservation of angular momentum (Equation (7)) and an inner radius R_X determined by magnetic truncation (Equation (36)). Given the expected disk radius and mass accretion rate onto the disk, the timescale for viscous evolution is shorter than that for disk mass accumulation provided that the dimensionless parameter α ≳ 10⁻⁵. As a result, the disk is expected to reach a steady state, with low mass compared to the central planet, and with a surface density distribution Σ(r) given by Equations (30) and (32).

Magnetic fields play an important role in determining the geometry of the inner disk and inner infall region. The balance between magnetic pressure and mass accretion rate through the disk defines the truncation radius R_X (Equation (36)), which specifies the inner disk edge. The field configuration also defines a capture radius R_σ (Equation (45)), which is determined by the balance between the ram pressure of infalling gas and the magnetic pressure at the point where the field lines are horizontal (see Figure 3).

The forming planet generates power from several contributions (see Section 5.1). The relevant luminosity sources include shocks from direct infall onto the planetary surface, accretion from the inner disk edge onto the planet, viscous accretion through the circumplanetary disk, as well as both shocks and mixing from material falling onto the disk surface. The planet generates additional internal power through gravitational contraction. The total luminosity is comparable to the benchmark power scale given by Equation (47).

With the properties of the components specified, we have calculated radiative signatures of the planet-forming process. The central planet and disk produce an intrinsic spectral energy distribution that corresponds to surface temperatures in the range T = 100–1000 K (Figure 4). The infalling envelope attenuates the central source luminosity, which is reradiated at longer wavelengths, primarily in the infrared. The resulting spectral energy distributions ν L_ν are significantly broader than single-temperature blackbodies and show systematic variations with planet mass M_p, accretion rate ${\dot{M}}_{{\rm{p}}}$, magnetic field strength B_p, and semimajor axis a (these trends are illustrated in Figures 5–9).

The infall-collapse-disk solutions constructed here for the formation of giant planets are analogous to those found earlier for the star formation problem (see the review of Shu et al. 1987). It is useful to outline the ways in which the formation processes are similar, and different, for the two cases. In both settings, angular momentum plays a key role in determining the geometry of the incoming material. Both star formation and planet formation produce nearly pressure-free collapse flow, with high-angular-momentum orbits that lead to the formation of a disk structure. Significantly, most of the mass initially falls onto the disk in both cases. The resulting dynamic range, as determined by the ratio of the disk radius to that of the central body, is much larger for star formation (R_C /R_* ∼ 10⁴) than for planet formation (R_C /R_p ∼100). Stars form much faster, with typical formation times of order 0.1 Myr, compared to 1 Myr for planets. ¹⁰ With longer formation times and smaller dynamic ranges, circumplanetary disks require much smaller viscosity to reach a steady state where accretion through the disk keeps pace with infall onto the disk. Although the source of viscosity remains under study, modest values of the parameter α ∼ 10⁻⁴ are sufficient to keep circumplanetary disks in steady state (and such values are expected to be realized). Much larger viscosity levels are required for star-forming disks to reach steady state. As a result, circumstellar disks are expected to experience gravitational instability during their early formative stages. In addition, the circumstellar disks are subject to episodic accretion, which produces FU Orionis outbursts. It remains to be seen if circumplanetary disks are subject to the same episodic behavior.

The dynamic range in temperature is also much larger for star formation, where the stellar surface has typical temperature T_* ∼ 6000 K, and the outer part of the infall region has interstellar temperatures T_ism ∼ 10−30 K. For planet formation, the planetary surface temperature is significantly smaller T_p ∼ 1000–1500 K, and the temperature at the outer boundary is typically T_H ∼ 100 K (for a ∼ 5 au). The temperature range for star formation (factor of ∼300) is more than an order of magnitude greater than the range for planet formation (factor of ∼12). Planets forming at smaller semimajor axes (e.g., a ∼ 0.1 au), have even smaller temperature ranges.

In both star and planet formation, the rate $\dot{M}$ at which the growing body gains mass is of fundamental importance. In the case of star formation, ${\dot{M}}_{\mathrm{core}}\approx {v}_{{\rm{s}}}^{3}/G$ is determined by the pre-collapse conditions in the background molecular cloud. This quantity determines the rate at which the entire star–disk system gains mass from its parental cloud. Dissipative processes within the circumstellar disk then determine the smaller rate ${\dot{M}}_{* {\rm{d}}}$ at which the disk accretes material onto the star. In the context of planet formation, the background circumstellar disk plays a role similar to the cloud in star formation. In the latter case, however, the semimajor axis of the planet determines the Hill radius, which in turn affects the accretion rate ${\dot{M}}_{{\rm{p}}}$ into the sphere of influence (where ${\dot{M}}_{{\rm{p}}}$ is a fraction of ${\dot{M}}_{* {\rm{d}}}$). Unlike stars, forming planets are thus subject to a systematic variation in background conditions, with corresponding variations in $\dot{M}$.

To help organize our understanding of the late stages of the planet formation process, we can consider the ordering of the relevant length scales:

$$\begin{eqnarray}&&{R}_{{\rm{p}}}\lt {R}_{X}\sim {R}_{\sigma }\ll {R}_{C}\sim {R}_{{\rm{d}}}\lt {R}_{H}\sim H\ll a.\end{eqnarray} \tag{ 72 }$$

During the phase when the growing planet accretes most of its gas, the protoplanet thus develops a hierarchical structure. Most of the mass entering the sphere of influence from the background circumstellar disk falls onto the circumplanetary disk, which is much larger than the planet itself. Magnetic fields lead to a truncation of the inner disk and provide a larger effective-capture cross section at a scale of several planetary radii. The nominal disk radius is one-third of the Hill radius, which is roughly comparable to the scale height of the circumstellar disk. For the formation of giant planets near or beyond the ice line (a ≳ 3 au), all of these scales are much smaller than the semimajor axis at the location of the planet. For hot Jupiters (a ≲ 0.1 au), however, the length scales that define the outer disk (R_d, R_C , R_H , H) become comparable to the magnetic radii (R_X , R_σ ). This reordering of length scales implies that giant planet formation could proceed differently in the inner and outer (circumstellar) disk, with a boundary between the two regimes at roughly a ∼ 0.1 au.

The relevant timescales also have a well-defined ordering:

$$\begin{eqnarray}&&{t}_{\mathrm{disk}}\lt {t}_{\mathrm{orb}}\ll {t}_{\mathrm{vis}}\ll {t}_{\mathrm{form}}\lesssim {t}_{\mathrm{kh}}.\end{eqnarray} \tag{ 73 }$$

The orbital timescale of the forming planet (around the star) is 9 times longer than the orbital timescale at the outer edge of the circumplanetary disk (around the planet), so that t_orb = 9 t_disk. The viscous evolution timescale t_vis = ${R}_{C}^{2}/\nu$ of the circumplanetary disk is much longer than the orbit time, but much shorter than the formation timescale t_form = ${M}_{{\rm{p}}}/{\dot{M}}_{{\rm{p}}}$. As a result, the disk can evolve to a steady-state configuration (see Equation (30)). The Kelvin–Helmholtz timescale t_kh sets the timescale for the evolution of the internal planetary structure, and is comparable to the formation time if the luminosity is determined by the scale L₀ from Equation (47), but somewhat longer if we use only the internally generated luminosity to define t_kh.

In general, the Hill radius R_H provides the outer boundary for the planetary sphere of influence, including in this present work. Although this ansatz provides a good first approximation, the picture is more nuanced. Numerical simulations indicate that some of the material that initially flows into the Hill sphere promptly flows back out (e.g., Lambrechts et al. 2019). The treatment of this paper accounts for this loss by taking the mass infall rate ${\dot{M}}_{{\rm{p}}}$ to include only the material that initially stays within the Hill sphere and falls toward the planet. A related complication is that the Hill radius is only an approximation to the true boundary. Pressure and magnetic corrections lead to smaller values for the effective boundary between the planet and its background disk (Appendix C). A smaller boundary radius, in turn, leads to a smaller centrifugal barrier R_C , which sets the nominal radius of the circumplanetary disk (Equations (5)–(7)). Two-dimensional simulations find that the expected disk radius at R_C = R_H /3 provides a good approximation (Martin & Lubow 2011), whereas three-dimensional simulations find somewhat smaller values for R_C (Fung et al. 2019). Three-dimensional simulations (e.g., Lambrechts et al. 2019) show that infall takes place preferentially along the poles of the system. This confinement results in a smaller centrifugal radius (${R}_{C}\to {R}_{C}[1-{\mu }_{b}^{2}]$) for a given planetary mass and a somewhat steeper surface density distribution. Finally, some of the material that initially falls onto the circumplanetary disk must eventually flow back out into the circumstellar disk. Viscous torques cause the circumplanetary disk to spread, so that most of the mass flows inward but some mass flows outward to conserve angular momentum. As a result, some fraction of the disk material must diffuse beyond the outer boundary (given by the generalized Hill radius) and rejoin the background circumstellar disk. ¹¹ The analytic treatment developed in this paper is robust, in that it provides solutions for any value of the centrifugal radius. On the other hand, additional work is required to determine the expected values of R_C to higher accuracy.

This analytic approach highlights how planet formation can vary with the location of the planetary orbit. For a given planetary mass, the physical size of the Hill radius and the circumplanetary disk radius grow with increasing semimajor axis a. Disk viscosity must play an increasingly important role for more distant planets. This trend has the effect of making planet formation (through the core-accretion paradigm) more difficult with increasing a. This trend acts in addition to the previously discussed effect that the formation of planetary cores takes longer for larger orbits. In the opposite limit, for small orbits akin to those of hot Jupiters, accretion through circumplanetary disks is even more efficient. However, for a ≲ 0.1 au, the expected magnetic truncation radius becomes smaller than the disk radius, and infall onto the entire planet–disk system can be suppressed.

This paper also provides an analytic description of the radiative signature of forming gas giants. Observations of dust emission from circumplanetary disks (which are embedded within circumstellar disks) have just become technologically possible (Benisty et al. 2021) and more such detections are expected. Current observations of circumplanetary disks are roughly consistent with theoretical expectations—including those of this paper. As future observations further elucidate the properties of forming planets, the analytic treatment of this paper can be readily adopted and generalized to provide theoretical descriptions of the formation process.